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Abstract 

Phenomenological plastic flow stress models are used extensively in the simulation of large deformations of 
metals at high strain-rates and high temperatures. Several such models exist and it is difficult to determine the 
applicability of any single model to the particular problem at hand. Ideally, the models are based on the underlying 
(subgrid) physics and therefore do not need to be recalibrated for every regime of application. In this work we 
compare the Johnson-Cook, Steinberg-Cochran-Guinan-Lund, Zerilli-Armstrong, Mechanical Threshold Stress, and 
Preston-Tonks-Wallace plasticity models. We use OFHC copper as the comparison material because it is well 
characterized. First, we determine parameters for the specific heat model, the equation of state, shear modulus 
models, and melt temperature models. These models are evaluated and their range of applicability is identified. We 
then compare the flow stresses predicted by the five flow stress models with experimental data for annealed OFHC 
copper and quantify modeling errors. Next, Taylor impact tests are simulated, comparison metrics are identified, and 
the flow stress models are evaluated on the basis of these metrics. The material point method is used for these 
computations. We observe that the all the models are quite accurate at low temperatures and any of these models 
could be used in simulations. However, at high temperatures and under high-strain-rate conditions, their accuracy 
can vary significantly. 



1 Introduction 

The Uintah computational framework Ide St. Germain et alJ (EoOO)) was developed to provide tools for the 
simulation of multi-physics problems such as the interaction of fires with containers, the explosive deformation and 
fragmentation of metal containers, impact and penetration of materials, dynamics deformation of air-filled metallic 
foams, and other such situations. Most of these situations involve high strain-rates. In some cases there is the 
additional complication of high temperatures. This work arose out of the need to validate the Uintah code and to 
quantify modeling errors in the subgrid scale physics models. 

Plastic flow stress models and the associated specific heat, shear modulus, melting temperature, and equation of state 
models are subgrid scale models of complex deformation phenomena. It is unreasonable to expect that any one model 
will be able to capture all the subgrid scale physics under all possible conditions. We therefore evaluate a number of 
models which are best suited to the regime of interest to us. This regime consists of strain-rates between 10 3 /s and 
10 6 /s and temperatures between 230 K and 800 K. We have observed that the comb ined effect of high temperatur e 
and high strain-rates has been g lossed over in most other similar works (for example Zeri Hi and Armstron 
Ijohnson and Holmauisll Jl988l) : IZocher et all J2000l) '). Hence we examine the temperature dependence of plastic 
deformation at high strain-rates in some detail in this paper. 
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In this paper, we attempt to quantify the modeling errors that we get when we model large-deformation plasticity (at 
high st rain-rates and high temper atures) with five recently developed models. T hese models are t he Joh nson-Cook 
model (Johnson and Cook ( 1983)), the Steinberg-Cochran-Guinan-Lund model (Steinberg et al. ( 1980); 
IStein berg and Lun dlJl989|)). the Zerilli- Arm strong model JZerilli and Armstrong! i 1987|)). the Mechanica l Threshold 
Stress model JFollansbee and KockH fl988a). an d the Pre ston-Tonk s-Wallace model dPreston et alJlEooH)). We a lso 
evalua te the associated shear modulus mod els oflVarshnil(ll970l).ISteinberg et al.Ul9 80l). andlNadal and Le Poad 
(2003). The melting temperature models of Steinberg et al. ( 1980) and Burakovskv etal.l f2000a) are also examined. 
A temperature-dependent specific heat relation is used to compute specific heats and a form of the Mie-Griineisen 
equation of state that assumes a linear slope for the Hugoniot curve are also evaluated. We suggest that the model that 
is most appropriate for a given set of conditions can be chosen with greater confidence once the modeling errors are 
quantified, 

The most common approach for determining modeling error is the comparison of predicted uniaxial stress-strain 
curves with experimental data. For high strain-rate conditions, flyer plate impact tests pr ovide further 
one-dimensional data that can be used to evaluate plasticity models. Taylor impact tests dTavloj lfl948)) can be use to 
obtain two-dimensional estimates of modeling errors. We restrict ourselves to comparing uniaxial tests and Taylor 
impact tests in this paper; primarily because high-temperature flyer plate impact experimental data are not readily 
available in the literature. We simulate uniaxial tests and Taylor impact tests with the Material Point Method 
dSulskvetal.Nl994ll995l) ). The model parameters that we use in these simulations are, for the most part, the values 
that are available in the literature. We do not recalibrate the models to fit the experimental data that we use for our 
comparisons. For simplicity, we use annealed OFHC copper as the material for which we evaluate all the models 
because this material is well -characterized. A similar exercise for various tempers of 4340 steel can be found 
elsewhere (Baneriee (2005a)). 

Most comparisons between experimental data and simulations involve the visual estimation of errors. For example, 
two stress-strain curves or two Taylor specimen profiles are overlaid on a graph and the viewer estimates the 
difference between the two. We extend this approach by providing quantitative estimates of the error and providing 
metrics with which such estimates can be made. The metrics are discussed and the models are evaluated on the basis 
of these metrics. 

The organization of this paper is as follows. Section[2]discusses the specific heat model, the equation of state, the 
melting temperature models, and the shear modulus models. Flow stress models are discussed in Section|5]and 
evaluated on the basis of one-dimensional tension and shear tests. Section|4]discusses experimental data, metrics, and 
simulations of Taylor impact tests. Conclusions are presented in Section|5] 



2 Models 

In most computations involving plastic deformation, the specific heat, the shear modulus, and the melting 
temperature are assumed to be constant. However, the shear modulus is known to vary with temperature and pressure. 
The melting temperature can increase dramatically at the large pressures experienced during high strain-rate 
deformation. In some materials, the specific heat can also change significantly with change in temperature. If the 
range of temperatures and strain-rates is small then these variations can be ignored. However, if a simulation involves 
a change in strain-rate from quasistatic to explosive, and a change in temperature from ambient values to values that 
are close to the melt temperature, the temperature- and pressure-dependence of these physical properties has to be 
taken into consideration. 

The models used in our simulations are discussed in this section. The material response is assumed to be isotropic. 
The stress is decomposed into a volumetric and a deviatoric part. The volumetric part of the stress is computed using 
the equation of state. The deviatoric part of the stress is computed using an additive decomposition of the rate of 
deformation into elastic and plastic parts, the von Mises yield condition, and a flow stress model. The variable shear 
modulus is used to update both the elastic and plastic parts of the stress and is also used by some of the flow stress 
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models. The melting temperature model is used to determine if the material has melted locally and also feeds into 
one of the shear modulus models. The increase in temperature due to the dissipation of plastic work is computed 
using the variable specific heat model. We stress physically -based models in this work because these can usually be 
used in a larger range of conditions than empirical models and need less recalibration. 

Copper shows significant strain hardening, strain-rate sensitivity, and temperature dependence of plastic flow 
behavior. The material is quite well characterized and a significant amount of experimental data are available for 
copper in the open literature. Hence it is invaluable for testing the accuracy of plasticity models and validating codes 
that simulate plasticity. In this work, we have only considered fully annealed oxygen-free high conductivity (OFHC) 
copper and electrolytic tough pitch (ETP) copper. 



2.1 Adiabatic Heating, Specific Heat, Thermal Conductivity 



A part of the plastic work done is converted into heat and used to update the temperature of a particle. The increase in 
temperature (AT) due to an increment in plastic strain ( Ae p ) is given by the equation 

AT = ^Ae p (1) 

where \ is the Taylor-Quinney coefficient ( Tavlor and Ou innevl CHI), a nd Cn is the specific heat. T he value of the 
Taylor-Quinney coefficient is assumed to be 0.9 in all our simulations (see lRavichandran et alJ J20011 for more 
details on the variation of \ with strain and strain-rate). The specific heat is also used in the estimation of the change 
in internal energy required by the Mie-Griineisen equation of state. 

The specific heat (C p ) versus temperature (T) model used in our simulations of copper has the form shown below. 
The units of C v are J/kg-K and the units of T are degrees K. 

_ fo. 0000416 T 3 — 0.027 T 2 + 6.21 T — 142.6 for T < 270K 
P ~ ] 0.1009 T + 358.4 for T > 270K 



A constant specific heat (usually assumed to be 414 J/kg-K) is not appropriate at temperatures below 250 K and 
temperatures above 700 K, as can be seen from Figure^ The specific heat predicted by our model (equation (|2ji) is 
shown as a solid line in the figure. This model is used to compute the specific heat in all the simulations described in 
this paper. 

The heat generated at a material point is conducted away at the end of a time step using the transient heat equation. 
The thermal conductivity of the material is assumed to be constant in our calculations. The effect of conduction on 
material point temperature is negligible for the high strain-rate problems simulated in this work. We have assumed a 
constant thermal conductivity of 386 W/(m-K) for copper which is the value at 500 K and atmospheric pressure. 



2.2 Equation of State 



The hydrostati c pressure (p ) is calculated us ing a tem perature-corrected Mie-Griineisen equation of state of the form 
used bv lZocher etail d200(I> (see also lWilkinsllll999h . p.61) 

w cg(,-i)[,-i»( r i)] , 
[x-s.(v-i)] 2 » 

where Co is the bulk speed of sound, p is the initial density, p is the current density, T is the Griineisen's gamma at 
reference state, S a — dU s /dU p is a linear Hugoniot slope coefficient, U s is the shock wave velocity, U p is the 
particle velocity, and E is the internal energy per unit reference specific volume. 
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The change in internal energy is computed using 



E=± [ Cv dT* Cv{T - To) (4) 
vo J Vq 

where Vq = 1/po is the reference specific volume at temperature T = T , and C v is the specific heat at constant 
volume. In our simulations, we assume that C p and C v are equal. 

The hydrostatic pressure is used to compute the volumetric part of the Cauchy stress tensor in our simulations. The 
parameters that we use in the Mie-Griineisen equation of state are shown in Tabled 

Figure[2]shows plots of the pressure predicted by the Mie-Griineisen equation of state at three different temperatures. 
The reference temperature for these calculations is 300 K. An initial density po of 8930 kg/m 3 has been used in the 
model calculations. The predicted pressures can be compared with pressures obtained from experimental shock 
Hugoniot data (shown by symbols in Figure^. The model equation of state performs well for compressions less than 
1.3. The pressures are underestimated at higher compression. We rarely reach compressions greater than 1.2 in our 
simulations. Therefore, the model that we have used is acceptable for our purposes. 



2.3 Melting Temperature 

The melting temperature model is used to determine the pressure-dependent melt temperature of copper. This melt 
temperature is used to compute the shear modulus and to flag the state (solid or liquid) of a particle. Two melting 
temperature models are evaluated in this paper. These are the Steinberg-Cochran-Guinan (SCG) melt model and the 
Burakovsky-Preston-Silbar (BPS) melt model. 



2.3.1 The Steinberg-Cochran-Guinan (SCG) melt model 

The Steinberg-Cochran-Guinan (SCG) melt model JSteinberg et alJ |l980)) is a relation between the melting 
temperature (T m ) and the applied pressure. This model is based on a modified Lindemann law and has the form 

arro-«-i/3) = £_ (5) 
Pa 

where T m o is the melt temperature at rj — 1, a is the coefficient of the first order volume correction to Griineisen's 
gamma (T ). 



T m (p) = T m0 exp 



2a 1 



2.3.2 The Burakovsky-Preston-Silbar (BPS) melt model 

An alternative melting relation that is based on dislo cation-m ediated phase transitions is the 
Burakovsky-Preston-Silbar (BPS) model ( B urakovskv et al j J200 0a)). The BPS model has the form 



T m (p) = T m (0) 



1 Po 



7i 7i 



4/3 



Po 



T m (0) 



nXpo v ws 



In 



87rln0- 1) k b \4& 2 j o c (7 1 , 



(6) 
(7) 



where p is the pressure, 77 is the compression (determined using the Murnaghan equation of state), /io is the shear 
modulus at room temperature and zero pressure, /i := dp /dp is the derivative of the shear modulus at zero pressure, 
Kq is the bulk modulus at room temperature and zero pressure, K := dK/ dp is the derivative of the bulk modulus 
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Figure 1 : Variation of the specific heat of coppe r with temperature. The s o lid line shows the values predicted b y the 
mode l. Symbols show experimen tal data from lOsborne an d Kirbv ( 1977), MacDonald and MacDonald (1981), and 
iDobrosavlievic and Madid i ll 99 ll) . 

Table 1: Parameters used in the Mie-Griin eisen EOS for copper. Th e bulk speed of sound and the slope of the 
linear fit to the Hugoniot for copper are from Mitchell and Nellis ( 1981). The value of the Griineisen gamma is from 
iMacDonald and MacDonaldl(ri98Tl). 

~c (m/s) s a r (r< 700K) r (r> 700K) 

3933 1.5 1.99 2.12 



1000 




Figure 2: The pressure predicted by the Mie-Griineisen equation of state for copper as a function of compression. The 
continuous lin es show the values predicted by the model for three temperatures. The symbols show exp erimental data 
obtained from Mc Queen et al J d 1 970l) . iMarshl ll 1 9 80l) . iMitchell and NelrisNl98ll) . and lWang et aT]d2000h . The original 
sources of the experimental data can be found in the above citations. 
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at zero pressure, k is a constant, A := b 3 /vws where b is the magnitude of the Burgers vector, vws is the 
Wigner-Seitz volume, z is the coordination number, a is a constant, p c (T m ) is the critical density of dislocations, and 
kb is the Boltzmann constant. 



2.3.3 Evaluation of melting temperature models 

Table|2]shows the parameters used in the melting temperature mod els of copper. F igure|3lshows a comparison of the 
two melting temperature models along with experimental data from lBurakovskv et al J ( 120 00a) (shown as open 
circles). An initial density po of 8930 kg/m 3 has been used in the model calculations. 

Both models predict the melting temperature quite accurately for pressures below 50 GPa. The SCG model predicts 
melting temperatures that are closer to experimental values at higher pressures. However, the data at those pressures 
are sparse and should probably be augmented before conclusions regarding the models can be made. In any case, the 
pressures observed in our computations are usually less than 100 GPa and hence either model would suffice. We have 
chosen to use the SCG model for our copper simulations because the model is more computationally efficient than 
the BPS model. 



2.4 Shear Modulus 

The shear modulus of copper decreases with temperature and is also pressure-dependent. The value of the shear 
modulus at room temperature is around 150% of the value close to melting. Hence, if we use the room temperature 
value of shear modulus for high temperature simulations we will overestimate the shear stiffness. This leads to the 
inaccurate estimation of the plastic strain-rate in radial return algorithms for elastic-plastic simulations. On the other 
hand, if the pressure-dependence of the shear modulus is neglected, modeling errors can accumulate for simulations 
involving shocks. 

Three model s for the shear m odulus (p) have been used in our simulations. The MTS shear modulus model was 
developed by Varshni|(|l970l) and ha s been used in conju nction with the Mechanical Threshold Stress (MTS) flow 
stress model JChen and GravMl996l):lGoto et alJfcOOOal) '). The Steinberg-Cochran-Guinan (SCG) shear modulus 
model was developed by Guinan and Steinberg ( 1974) and has been used in conjunction with the 
Steinberg -Cochran-Guinan -Lund (SCGL) flow stress model. The Nadal and LePoac (NP) shear modulus model 
(Nadal and Le Poac is a recently developed model that uses Lindemann theory to determine the temperature 

dependence of shear modulus and the SCG model for pressure dependence. 



2.4.1 MTS Shear Modulus Model 



The MTS shear modulus model has the form JVarshnilfl97(il) : IChen and Gravl fl99 



H(T) = po — — - (8) 

exp{T /T) - 1 

where po is the shear modulus at OK, and D, To are material constants. The shortcoming of this model is that it does 
not include any pressure-dependence of the shear modulus and is probably not applicable for high pressure 
applications. However, the MTS shear modulus model does capture the flattening of the shear modulus-temperature 
curve at low temperatures that is observed in experiments. 
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Table 2 : Parameters used in melti ng temperature models for copper. The parameter T m n used in the SCG model 
is from lGuinan and Steinberg! Jl974l) . The value of To is from lMacDonald and MacDonaldl Jl98ll) . The value of a 
has been chosen to fit the experiment al data. The values of the initial bulk and shear moduli and their derivatives 
in the BPS model are fronuG uinan and Steinberg ( 1974). The remaining parameters for the BPS model are from 
iBurakovskv and Prestonl J2000h and lBurakovskv et alJ J2000bl) . 

Steinberg-Cochran-Guinan (SCG) model 

T m0 (K) Tp a 

1356.5 1.99 1.5 

Burakovsky-Preston-Silbar (BPS) model 

Ko (GPa) K (GPa) fj, k z b 2 p c (T m ) a A vws a (nm) 
~ 137 5.48 47.7 1.4 1.25 12 0.64 2.9 1.41 a 3 /4 3.6147 
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-50 50 100 150 200 

Pressure (GPa) 

Figure 3: The melting temperature of copper as a function of pressure. The lines sho w values predicted by the SCG 
and BPS models. The open circles show experimental data obtained from Bura kovskv et alJ Il2000ah . The original 
sources of the experimental data can be found in the above citation. 



O Burakovsky et al. (2000) 

SCG Melt Model 

BPS Melt Model 
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2.4.2 SCG Shear Modulus Model 



The Steinberg-Cochran-Guinan (SCG) shear modulus model JSteinberg et alJ Jl980l) : lzocher et all J2000l) 1 is 
pressure dependent and has the form 

M ( P) T) = Mo + ^^3+^(T-300); r, = P / Po (9) 

where, (iq is the shear modulus at the reference state(T = 300 K, p = 0, r\ = 1), p is the pressure, and T is the 
temperature. When the temperature is above T m , the shear modulus is instantaneously set to zero in this model. 



2.4.3 NP Shear Modulus Model 



The Nadal-Le Poac (NP) shear modulus model ( Nadal and Le Poad i200$> ) is a modified version of the SCG model. 
The empirical temperature dependence of the shear modulus in the SCG model is replaced with an equation based on 
Lindemann melting theory. In addition, the instantaneous drop in the shear modulus at melt is avoided in this model. 
The NP shear modulus model has the form 



i 



J(T) 



Mo 



d\x p 
dp 77 1 / 3 



(1 — T) + — f— k b T 

Cm 



C := 



(10) 



where 



J(f) := 1 + exp 



1 + 1/C 



1 + C/(1-T) 



for f : =JL e [0,l + C], 

J- m. 



(ID 



/z is the shear modulus at K and ambient pressure, £ is a material parameter, is the Boltzmann constant, m is the 
atomic mass, and / is the Lindemann constant. 



2.4.4 Evaluation of shear modulus models 

The parameters used in the three shear modulus models are given in Table[5] Figure|4la) shows the shear modulus 
predicted by the MTS shear modulus model at zero hydrostatic pressure. It can be seen that the model fits the low 
temperature data quite well. The shear moduli predicted by the SCG and NP shear models are shown in Figure|4jb) 
and Figure|3c), respectively. The SCG shear model predicts slightly different moduli than the NP model at different 
values of compression. Both models fit the experimental data quite well except at very low temperatures (at which the 
MTS model performs best). We have not be able to validate the pressure dependence of the shear modulus at high 
temperatures due to lack of experimental data. An initial density of 8930 kg/m 3 has been used in the model 
calculations. 



3 Flow Stress Models 

We have explored five temperature and strain-rate dependent models that can be used to compute the flow stress: 

1 . the Johnson-Cook model 

2. the Steinberg-Cochran-Guinan-Lund model. 

3. the Zerilli-Armstrong model. 
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Table 3: Parameters used in shear modulus models for copper. The parameter s for the M TS model have been chosen 
to fit the experiment al data. The parame t ers for the SCG model are from lGuinan and SteinberdJl974l) . The NP model 
parameters are from lNadal and Le PoacH2003l) . 



MTS shear modulus model 




SCG shear modulus model 


Mo (GPa) D (GPa) T (K) 




Mo (GPa) dfi/dp dp/dT(GPa/K) 


51.3 3.0 165 




47.7 1.3356 0.018126 


NP shear modulus model 


Mo (GPa) dfi/dp C 


C 


m (amu) 


50.7 1.3356 0.04 


0.057 


63.55 




n Overton et al. (1955) 

O Nadal and LePoac (2003) 




0.2 0.4 0.6 0.8 1 

T/T 
m 



(c) NP Shear Model 



Figure 4: Shear modulu s of copper as a f unction of temperature and pressure. The symbols represent experimental 
data from lOverton and GaffnevMl955l) and lNadal and Le Poad J2003I) . The lines show values of the shear modulus at 
different compressions (rj — p/po). 
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4. the Mechanical Threshold Stress model. 

5. the Preston-Tonks-Wallace model. 



The Johnson-Cook (JC) model Jjohnson and Cookl i 19831) 1 is purely empirical and is the most widely used of the five. 
However, this model exhibits an unrealistically s mall strain-rate dependence at high temperatures. Th e 
Steinberg-Cochran-Guinan-Lund (SCGL) model JSteinberg et al.Nl98(]l) : ISteinberg and Lundl i ll 9891V ) is 
semi-empirical. Th e model is purely empirical a nd strain-rate independent at high strain-rates. A dislocation-based 
extension based on Hoge and Mukherie e Nl977l) is used at low stra in-rates. The SCGL model is used extensively by 
the shock physics community. The Zerilli- Armstrong (ZA) model (Zerilli and Armstrong lll987t) 1 is a simple 
physically-based model that has been used extensively. A more complex mo del that is based on ideas from 
dislocation dynamics is the Mechanical Threshold Stress (MTS) model jFollansbee and Kocksl £l988)). This model 
has been used to model the plastic def ormation of copper, tant alum JChen an TGrayfT l996l) '). alloys of steel 
jGoto et alJ j2000al) : lBanerieel(l200 5a)). and aluminum alloys ( Puc hi-Cabrera et alJ J2QQ lb S . Howe ver, the MTS 
model is limited to strain-rates less than around 10 7 /s. The Preston-Tonks-Wallace (PTW) model JPreston et alJ 
(120031) 1 is also physically based and has a form similar to the MTS model. However, the PTW model has components 
that can model plastic deformation in the overdriven shock regime (strain-rates greater that 10 7 /s). Hence this model 
is valid for the largest range of strain-rates among the five flow stress models. 



3.1 JC Flow Stress Model 

The Johnson-Cook (JC) model (Johns on and Cookl dl983l) 1 is purely empirical and gives the following relation for the 
flow stress (a y ) 

a y (e p ,e p ,T) = [A + B(e p y l ][l + CHe;)}[l-(T*) m ] (12) 
where e p is the equivalent plastic strain, e p is the plastic strain-rate, and A, B, C, n, m are material constants. 

The normalized strain-rate and temperature in equation fl!2l > are defined as 

e» (T - T ) 
? p :=JL an d T* := j- £L (13) 

where e p q is a user defined plastic strain-rate, Tq is a reference temperature, and T m is a reference melt temperature. 
For conditions where T* < 0, we assume that m = 1. 



3.2 SCGL Flow Stress Model 

The Stein berg-Cochra n-Guinan-Lund (SCGL) model is a semi-empirical model that was developed by 
Steinberg et al Idl980l) for high strain-rate situations and extended to low strain-rates and bcc materials by 
sTdribergari^Lundl ( ll989l) . The flow stress in this model is given by 

j^j /. (T a f and at < a p (14) 

Mo 

where a a is the athermal component of the flow stress, f(e p ) is a function that represents strain hardening, tr t is the 
thermally activated component of the flow stress, fx(p, T) is the pressure- and temperature-dependent shear modulus, 
and Mo is the shear modulus at standard temperature and pressure. The saturation value of the athermal stress is <7 ma x- 
The saturation of the thermally activated stress is the Peierls stress (a p ). The shear modulus for this model is usually 
computed with the SCG shear modulus model. 

The strain hardening function (/) has the form 

f(e p ) = [l + 0(e p + e pi )] n (15) 
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where (3, n are work hardening parameters, and e P i is the initial equivalent plastic strain. 

The thermal component (at) is computed using a bisection algorithm from the following equation 
(citetHoge77,Steinberg89). 



^ 6XP 



2U k 
k h T 



C2 
0* 



at < a v 



(16) 



where 2Uk is the energy to form a kink-pair in a dislocation segment of length Ld, kb is the Boltzmann constant, a p is 
the Peierls stress. The constants G\, C2 are given by the relations 



C x := 



p d L d ab 2 v 



Co 



D 
~P~6& 



(17) 



where pd is the dislocation density, L d is the length of a dislocation segment, a is the distance between Peierls 
valleys, b is the magnitude of the Burgers' vector, v is the Debye frequency, w is the width of a kink loop, and D is 
the drag coefficient. 



3.3 ZA Flow Stress Model 

The Zerilli-Armstrong (ZA) model JZerilli and Armstrong! ( fl987l 1 1 993h : Izerilljl 12004 )) is based on simplified 
dislocation mechanics. The general form of the equation for the flow stress is 

a v {e p ,e p ,T) = a a + B exp(-/3(e p )T) + B ^/e^cxp(-a(e p )T) . (18) 

In this model, a a is the athermal component of the flow stress given by 

a a :=a g + ^ + Ke$, (19) 

where a g is the contribution due to solutes and initial dislocation density, kh is the microstructural stress intensity, I is 
the average grain diameter, K is zero for fee materials, B, Bo are material constants. 

In the thermally activated terms, the functional forms of the exponents a and (3 are 

a = a - 011 ln(ep); /3 = fa - /3i ln(e p ); (20) 

where ao, 0:1, 0o, (3i are material parameters t hat depend on the type of m aterial (fee, bec, hep, alloys). The 
Zerilli-Armstrong model has been modified by Abed and Voviadiis (2005) for better performance at high 
temperatures. However, we have not used the modified equations in our computations. 



3.4 MTS Flow Stress Model 

The Mechanical Threshold Stress (MTS) model jFollansbee and Kocksl i ll 988b : lGoto et al.l ( l2000bl) : lKocksl(l200lh 
has the form 

a y (e p ,e p ,T) = a a + {S. l a. l + S e a e ) ^ P,T ^ (21) 

Mo 

where a a is the athermal component of mechanical threshold stress, cr, is the component of the flow stress due to 
intrinsic barriers to thermally activated dislocation motion and dislocation-dislocation interactions, a e is the 
component of the flow stress due to microstructural evolution with increasing deformation (strain hardening), (Si, S e ) 
are temperature and strain-rate dependent scaling factors, and /iq is the shear modulus at K and ambient pressure, 
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The scaling factors take the Arrhenius form 



5; 



k h T 



S e 



g 0i b 3 n{p, T) 

k b T 
g 0e b 3 fi(p, T) 



In 



In 



e pOc 



1/9, 



1/Pi 



1/Pe 



(22) 



(23) 



where fcf, is the Boltzmann constant, 6 is the magnitude of the Burgers' vector, (goi, <?0e) are normalized activation 
energies, (e p oi, e p o e ) are constant reference strain-rates, and (qi,p%, q e ,Pe) are constants. 



The strain hardening component of the mechanical threshold stress (er e ) is given by an empirical modified Voce law 

K<Te) (24) 



d<7 £ 
de„ 



where 



70 = 



9 Q [1-F(a e ) 
a + a\ In e p 



f 9 IV F(<j e ) 

Cl2\/%- CL3T 



ln(- 



O"0e 



tanh a — 

tanh(a) 
kT 

g 0es b 3 n(p,T) 



In 



(25) 
(26) 

(27) 
(28) 



and #0 is the hardening due to dislocation accumulation, 9jy is the contribution due to stage-IV hardening, 
(do, ax, »2, 03 1 a) are constants, <r es is the stress at zero strain hardening rate, <j 0es is the saturation threshold stress 
for deformation at K, go es is a constant, and e p oes is the maximum strain-rate. Note that the maximum strain-rate is 
usually limited to about 10 7 /s. 



3.5 PTW Flow Stress Model 



The Preston-Tonks-Wallace (PTW) model JPreston et alJ d2003l) ^ attempts to provide a model for the flow stress for 
extreme strain-rates (up to 10 n /s) and temperatures up to melt. A linear Voce hardening law is used in the model. 
The PTW flow stress is given by 



t s + a In 



2t 8 h(p,T) 




with 



a :- 



so 



1 — ip exp — j3 







H{p, T) thermal regime 
shock regime 



tp := exp(/3) - 1 



(29) 



(30) 



d a 

where r s is a normalized work-hardening saturation stress, s is the value of t s at OK, r y is a normalized yield stress 
9 is the hardening constant in the Voce hardening law, and d is a dimensionless material parameter that modifies the 
Voce hardening law. 
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The saturation stress and the yield stress are given by 



t s = max | s - (so - Soo)erf 
T y = max i y - (y - y^erf 



kT In 



it 



. s 



IJ-2 



P 



mm < yi — 7 , s — ■ 



7^, 



(31) 
(32) 



where Soc is the value of r s close to the melt temperature, (yo, j/oo) are the values of r y at OK and close to melt, 
respectively, (k, 7) are material constants, T — T/T m , (s 1; y 1; y 2 ) are material parameters for the high strain-rate 
regime, and 

/ \ 1 / 3 / \\ 1 / 2 

1 / Airp\ I (x(p, T) \ 



^ 2 I 3M 



P 



(33) 



where p is the density, and M is the atomic mass. 



3.6 Evaluation of flow stress models 

In this section, we evaluate the flow stress models on the basis of one-dimensi onal tension and compre ssion tests. The 
high rate tests have been simulated using the explicit Material Point Method Sulskv "et"al] l l 19941 1995b (see 
AppenedixIXl in conjunction with the stress update algorithm given in Append ix|5| The quasistatic tests have been 
simulated with a fully implicit version of the Material Point Method dGuilkev and Weissl J2003I) ') with an implicit 
stress update dSimo and HugheslJl998l) ). Heat conduction is performed at all strain-rates. As expected, we obtain 
nearly isothermal conditions for the quasistatic tests and nearly adiabatic conditions for the high strain-rate tests. We 
have used a constant thermal conductivity of 386 W/(m-K) for copper which is the value at 500 K and atmospheric 
pressure. To damp out large oscillatio ns in h igh strain-rate tests, we use a three-dimensional form of the von 
Neumann artificial viscosity ( Wilkins ( 1999), p. 29). The viscosity factor takes the form 

q = C oP l ] j^\trD\ + C 1 pl 2 (irDf (34) 

where Co and C\ are constants, p is the mass density, K is the bulk modulus, D is the rate of deformation tensor, and 
I is a characteristic length (usually the grid cell size). We have used Co = 0.2 and Cl = 2.0 in all our simulations. 
The temperature-dependent specific heat model, the Mie-Griineisen equation of state, and the SCG melting 
temperature model have been used in all the following simulations. 

The predicted s t ress -strain curves are compared with experimental dat a for a nnealed OFHC copper from tension tests 
dNemat-Nasserl J2004I) (p. 241-242)) and compression tests (Samantal (ll97ll) '). The data are presented in form of true 
stress versus true strain. Note that detailed verification has been performed to confirm the correct implementation of 
the models withing the Uintah code. Also note that the high strain-rate experimental data are suspect for strains less 
than 0.1. This is because the initial strain-rate fluctuates substantially in Kolsky-Hopkinson bar experiments. 



3.6.1 Johnson-Cook Model. 

The parameters that we have used in the Johnson-Cook (JC) flow stress model of annealed copper are given in 
Table |4] We have used the NP shear modulus model in simulations involving the JC model. 

The Johnson-Cook model is independent of pressure. Hence, the predicted yield stress is the same in compression 
and tension. The use of a variable specific heat model leads to a reduced yield stress at 77 K for high strain rates. 
However, the effect is relatively small. At high temperatures, the effect of the higher specific heat is to reduce the rate 
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of increase of temperature with increase in plastic strain. This effect is also small. The temperature dependence of the 
shear modulus does not affect the yield stress. However, it has a small effect on the value of the plastic strain-rate. 

The solid lines in Figures|5Ja) and (b) show predicted values of the yield stress for various strain-rates and 
temperatures. The symbols show the experimental data. The Johnson-Cook model overestimates the initial yield 
stress for the quasistatic (0. 1/s strain-rate), room temperature (296 K), test. The rate of hardening is underestimated 
by the model for the room temperature test at 8000/s. The strain-rate dependence of the yield stress is underestimated 
at high temperature (see the data at 1 173 K in Figure|3a)). For the tests at a strain-rate of 4000/s (Figure|3b)), the 
yield stress is consistently underestimated by the Johnson-Cook model. 

3.6.2 Steinberg-Cochran-Guinan-Lund Model. 

The parameters used in the Steinberg-Cochran-Guinan-Lund (SCGL) model of annealed OFHC copper are listed in 
Table|5] We have used the SCG shear modulus model in simulations involving the SCGL model. We could 
alternatively have used the NP shear modulus model. However, we use the SCG model to highlight a problem with 
the equivalence of dfi/dT and da y /dT that is assumed by the SCGL model. A bisection algorithm is used to 
determine the thermally activated part of the flow stress for low strain-rates (less than 1000/s). 

The solid lines in Figures[6ja) and (b) show the flow stresses predicted by the SCGL model. Clearly, the softening 
associated with increasing temperature is underestimated by the SCGL model though the yield stress at 8000/s is 
predicted reasonably accurately. For the tests at 4000/s shown in Figure|5Ib), the SCGL modes performs 
progressively worse with increasing temperature. 

Overall, at low temperatures, the high strain-rate predictio ns from the SC GL model match the experimental data best. 
This is not surprising since the original model by ISteinberg et ail ( ll980l) (SCG) was rate-i ndependent and designed 
for high strain-rate applications. However, the low strain rate extension by ISteinberg and Lundl ll989) does not lead 
to good predictions of the yield stress of OFHC copper at low temperatures. 

The high temperature response of the SCGL model is dominated by the shear modulus model; in particular, the 
derivative of the shear modulus with respect to temperature. From Figure fflb) we c an see that a value of -0.018126 
GPa/K for dfi/dT matches the experimental data quite well. Steinbe rg et alJ ([1980) assume that the values of 
(da y /dT)/<7yo and (dfi/dT) / fi Q (-3.8 x 10— 4 /K) are comparable. That does not appear to be the case for OFHC 
copper. 

If we extract the yield stresses at a strain of 0.2 from the experimental data shown in Figure[6jb), we get the following 
values of temperature and yield stress for a strain-rate of 4000/s: (77 K, 380 MPa); (496 K, 300 MPa); (696 K, 230 
MPa); (896 K, 180 MPa); (1096 K, 130 MPa). A straight line fit to the data shows that the value of da y /dT is -0.25 
MPa/K. The yield stress at 300 K can be calculated from the fit to be approximately 330 MPa. This gives a value of 
-7.6x 10~4 /K for (dcr y / dT) / a y o\ approximately double the slope of the shear modulus versus temperature curve. 
Hence, a shear modulus derived from a shear modulus model canno t be used as a m ultiplier to the yield stress in 
equation dl4> . Instead, the original form of the SCG model ( St einberg et alJ Jl980i) ) must be used, with the term 
(dfi/dT) /hq replaced by (da y /dT) /a y0 in the expression for yield stress. 

Figures^a) and (b) show the predicted yield stresses from the modified SCGL model. These plots show that there is 
a considerable improvement in the prediction of the temperature dependence of yield stress if the value of 
(da y /dT)/a y o is used instead of (dfi/dT)/ fiQ. However, the strain-rate dependence of OFHC copper continues to 
be poorly modeled by the SCGL model. 
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Table 4: Parameters used in the Johnson-Cook model for copper d.Tohnson and Cookl {T985)). 



A (MPa) 


B (MPa) 


C n 


m e p0 (/s) 


T (K) 


T m (K) 


90 


292 


0.025 0.31 


1.09 1.0 


294 


1356 



OFHC Copper (Johnson-Cook) 
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0.1 /s, 296K (Expt.) 
0.1 /s, 296K (Sim.) 
8000/s, 296K (Expt.) 
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2300/s, 873K (Sim.) 
1800/s, 1023K(Expt.) 
1800/s, 1023K(Slm.) 
0.066/s, 1173K(Expt.) 
0.066/s, 1173K(Sim.) 
960/s, 1173K(Expt.) 
960/s, 1173K(Slm.) 



0.4 0.6 
True Strain 

(a) Various strain-rates and temperatures. 
OFHC Copper (Johnson-Cook) 
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4000/s, 77K (Expt.) 
4000/s, 77K (Sim.) 

a 4000/s, 496K (Expt.) 
4000/s, 496K (Sim.) 

□ 4000/s, 696K (Expt.) 
4000/s, 696K (Sim.) 

4000/s, 896K (Expt.) 
4000/s, 896K (Sim.) 

v 4000/s, 1096K (Expt.) 
4000/s, 1096K (Sim.) 











0.2 



0.4 0.6 
True Strain 



0.8 



(b) Various temperatures at 4000/s strain-rate. 

Figure 5: Predict ed values of yiel d stress from the Johnson-Cook model. The experimental data at 873 K, 1023 K, and 
1173 K are from Samant aMl97ll) and represent compression tests. The remaining experimental data are from tension 
tests in Ne mat-Nassen (120 04 ) . The solid lines are the predicted values. 
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Table 5: Parameters used in the Steinberg-Cochran-Guinan-Lund model for copper. The parameters for the athermal 
part of the SCGL model are from ISteinberg et alJ dl98CA . The parameters fo r the thermally activated part of the model 
are from a number of sources. The estimate for the Peierls stress is based on lHobartl dl965h . 



a a (MPa) 


<x max (MPa) 


e pi (/s) n Ci (/s) 


U k (eV) 


a p (MPa) 


C* 2 (MPa-s) 
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8000/s, 296K (Expt.) 
8000/s, 296K (Sim.) 
1800/s, 1023K(Expt.) 
1800/s, 1023K(Sim.) 
0.066/s, 1173K(Expt.) 
0.066/s, 1173K(Sim.) 
960/s, 1173K(Expt.) 
960/s, 1173K(Sim.) 
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(a) Various strain-rates and temperatures. 
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4000/s, 77K (Sim.) 

a 4000/s, 496K (Expt.) 
4000/s, 496K (Sim.) 

□ 4000/s, 696K (Expt.) 
4000/s, 696K (Sim.) 

« 4000/s, 896K (Expt.) 
4000/s, 896K (Sim.) 

v 4000/s, 1096K (Expt.) 
4000/s, 1096K (Sim.) 
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0.4 0.6 
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(b) Various temperatures at 4000/s strain-rate. 



Figure 6: Predicted values of yield stress from the Steinberg-Cochran-Guinan-Lund model. Please see the caption of 
Figure|5]for the sources of the experimental data. 
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OFHC Copper (Steinberg-Cochran-Guinan-Lund) 
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4000/s, 696K (Sim.) 
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(b) Various temperatures at 4000/s strain-rate. 

Figure 7: Predicted values of yield stress from the modified Steinberg-Cochran-Guinan-Lund model. Please see the 
caption of Figure[5]for the sources of the experimental data. 
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3.6.3 Zerilli-Armstrong Model. 



In contrast to the Johnson-Cook and the Steinberg-Cochran-Guinan models, the Zerilli-Armstrong (ZA) model for 
yield stress is based on dislocation mechanics and hence has some physical basis. The parameters used for the ZA 
model are given in Tabled We have used the NP shear modulus model in our simulations that involve the ZA model. 

Figures [SJa) and (b) show the yield stresses predicted by the ZA model. From Figure^a), we can see that the ZA 
model predicts the quasistatic, room temperature yield stress quite accurately. However, the room temperature yield 
stress at 8000/s is underestimated. The initial yield stress is overestimated at high temperatures; as are the saturation 
stresses. 

Stress-strain curves at 4000/s are shown in Figure[8lb). In this case, the ZA model predicts reasonable initial yield 
stresses. However, the decrease in yield stress with increasing temperature is overestimated. We notice that the 
predicted yield stress at 496 K overlaps the experimental data for 696 K, while the predicted stress at 696 K overlaps 
the experimental data at 896 K. 

3.6.4 Mechanical Threshold Stress Model. 

The Mechanical Threshold Stress (MTS) model is different from the three previous models in that the internal 
variable that evolves in time is a stress (a e ). The value of the internal variable is calculated for each value of plastic 
strain by integrating equation J24i along a constant temperature and strain-rate path. An unconditionally stable and 
second-order accurate midpoint integration scheme has been used to determine the value of a e . Alternatively, an 
incremental update of the internal variable could be done using quantities from the previous timestep. The integration 
of the evolution equation is no longer along a constant temperature and strain-rate path in that case. We have found 
that two alternatives give us similar values of a e in the simulations that we have performed. The incremental update 
of the value of a e is considerably faster than the full update along a constant temperature and strain-rate path. 

The parameters for the MTS model are shown in Table0 The pressure-independent MTS shear modulus model has 
been used in simulations that use the MTS flow stress model. The reason for this choice is that the parameters of the 
model have been fit with such a shear modulus model. If the shear modulus model is changed, certain parameters of 
the model will have to be changed to reflect the difference. 

Figures|9ja) and (b) show the experimental values of yield stress for OFHC copper versus those computed with the 
MTS model. 

From Figure|9ja), we can see that the yield stress predicted by the MTS model almost exactly matches the 
experimental data at 296 K for a strain-rate of 0.1/s. The yield stress for the test conducted at 296 K and at 8000/s is 
underestimated. Though reasonably accurate yield stresses are predicted at 1023 K and 1800/s, the experimental 
curves exhibit earlier saturation than the model predicts. The same is true at 873 K and 2300 /s. The predicted yield 
stress is higher for the quasistatic test at 1 173 K than that observed experimentally. However, the higher rate test at the 
same temperature matches the experiments quite well except for a higher amount of strain hardening at large strains. 

The variation of yield stress with temperature at a strain-rate of 4000/s is shown in Figure|9[b). The figure shows that 
the yield stress is underestimated by the MTS model at all temperatures except 1096 K. The experimental data shows 
stage III or stage IV hardening which is not predicted by the MTS model that we have used. 

3.6.5 Preston-Tonks- Wallace Model. 

The Preston-Tonks-Wallace (PTW) model attempts to provide a single approach to model both thermally activated 
glide and overdriven shock regimes. The overdriven shock regime includes strain-rates greater than 10 7 . The PTW 
model, therefore, extends the possibility of modeling plasticity beyond the range of validity of the MTS model. We 
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Table 6: Parameters used in the Zerilli-Armstrong model for copper dZerilli and Armstrong! lll987l) ). 



u g (MPa) 


k h (MPa-mm 1 / 2 ) 
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0.1 /s, 296K (Expt.) 
0.1 /s, 296K (Sim.) 
8000/s, 296K (Expt.) 
8000/s, 296K (Sim.) 
2300/s, 873K (Expt.) 
2300/s, 873K (Sim.) 
1800/s, 1023K(Expt.) 
1800/s, 1023K(Sim.) 
0.066/s, 1173K(Expt.) 
0.066/s, 1173K(Sim.) 
960/s, 1173K(Expt.) 
960/s, 1173K(Sim.) 
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(a) Various strain-rates and temperatures. 
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□ 4000/s, 696K (Expt.) 
4000/s, 696K (Sim.) 

4000/s, 896K (Expt.) 
4000/s, 896K (Sim.) 

v 4000/s, 1096K (Expt.) 
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(b) Various temperatures and 4000/s strain-rate. 



Figure 8: Predicted values of yield stress from the Zerilli-Armstrong model. The symbols represent experimental data. 
The solid lines represented the computed stress-strain curves. Please see the caption of Figure[5]for the sources of the 
experimental data. 
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Table 7: Parameters used in the Mechanical Threshold Stress model for copper dFollansbee and Kocksl {l988)). 



° a (MPa) 


b (nm) 


<Ti (MPa) 


9m 


£p0i V s ) 


Pi 




40 


0.256 





1 


1 


1 


1 
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Pe 




er es (MPa) 
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ai (MPa-log(s)) 
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OFHC Copper (Mechanical Thresold Stress) 



700 




0.1 /s, 296K (Expt.) 
0.1 /s, 296K (Sim.) 
8000/s, 296K (Expt.) 
8000/s, 296K (Sim.) 
2300/s, 873K (Expt.) 
2300/s, 873K (Sim.) 
1800/s, 1023K(Expt.) 
1800/s, 1023K(Sim.) 
0.066/s, 1173K(Expt.) 
0.066/s, 1173K(Sim.) 
960/s, 1173K(Expt.) 
960/s, 1173K(Sim.) 
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(a) Various strain-rates and temperatures. 
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4000/s, 
4000/s, 
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4000/s, 
4000/s, 
4000/s, 
4000/s, 



77K (Expt.) 
77K (Sim.) 
496K (Expt.) 
496K (Sim.) 
696K (Expt.) 
696K (Sim.) 
896K (Expt.) 
896K (Sim.) 
1096K (Expt.) 
1096K (Sim.) 



0.4 0.6 
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(b) Various temperatures at 4000/s strain-rate. 



Figure 9: Predicted values of yield stress from the Mechanical Threshold Stress model. Please see the caption of 
Figure|5]for the sources of the experimental data. 
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have not conducted a simulations of overdriven shocks in this paper. However, the PTW model explicitly accounts for 
the rapid increase in yield stress at strain rates above 1000 /s. Hence the model is a good candidate for the range of 
strain-rates and temperatures of interest to us. The PTW model parameters used in our simulations are shown in 
Table|8] In addition, we use the NP shear modulus model in all simulations involving the PTW yield stress model. 

Experimental yield stresses are compared with those predicted by the PTW model in Figures [ToF a) and (b). The solid 
lines in the figures are the predicted values while the symbols represent experimental data. From Figure [Tol a) we can 
see that the predicted yield stress at 0.1/s and 296 K matches the experimental data quite well. The error in the 
predicted yield stress at 296 K and 8000/s is also smaller than that for the MTS flow stress model. The experimental 
data at 873 K, 1023 K, and 1 173 K were used by Preston et al. (2003) to fit the model parameters. Hence it is not 
surprising that the predicted yield stresses match the experimental data better than any other model. 

The temperature-dependent yield stresses at 4000/s are shown in Figure [TUf b). In this case, the predicted values at 77 
K are lower than the experimental values. However, for higher temperatures, the predicted values match the 
experimental data quite well for strains less than 0.4. At higher strains, the predicted yield stress saturates while the 
experimental data continues to show a significant amount of hardening. The PTW model predicts better values of 
yield stress for the compression tests while the MTS model performs better for the tension tests. 

3.6.6 Errors in the flow stress models. 

In this section, we use the difference between the predicted and the experimental values of the flow stress as a metric 
to compare the various flow stress models. The error in the true stress is calculated using 



A detailed discussion of the differ ences be tween the predicted and experimental true stress for one-dimensional tests 
can be found elsewhere (Baneriee (2005b)). In this paper we summarize these differences in the form of error 
statistics as shown in Tablesl9landllOI Only true strains greater than 0.1 have been considered in the generation of 
these statistics. The statistics in Tablesl9landll0lclearlv show that no single model is consistently better than the other 
models under all conditions. 

We can further simplify our evaluation by considering a single metric that encapsulates much of the information in 
these tables. Tablell llshows comparisons based on one such simplified error metric. We call this metric the average 
maximum absolute (MA) error. The maximum absolute (MA) error is defined as the sum of the absolute mean error 
and the standard deviation of the error. The extreme values of the error are therefore ignored by the metric and only 
values that are within one standard deviation of the mean are considered. 

From Tablell II we observe that the least average MA error for all the tests is 17% while the greatest average MA error 
is 64%. The PTW model performs best while the SCGL model performs worst. In order of increasing error, the 
models may be arranged as PTW, MTS, ZA, JC, and SCGL. 

If we consider only the tension tests, we see that the MTS model performs best with an average MA error of 14%. 
The Johnson-Cook model does the worst at 25% error. For the compression tests, the PTW model does best with an 
error of 10% compared to the next best, the MTS model with a 35% error. The SCGL error shows an average MA 
error of 126% for these tests. 

For the high strain-rate tests, the MTS model performs better than the PTW model with an average MA error of 15% 
(compared to 18% for PTW). The low strain-rate tests are predicted best by the PTW model (5 %) and worst by the 
SCGL model (219 %). Note that this average error is based on two tests at 296 K and 1 173 K and may not be 
representative for intermediate temperatures. 

The PTW model shows an average MA error of 16% for the high temperature tests compared to 27% for the MTS 




(35) 
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Table 8: Parameters used in the Preston-Tonks-Wallace yield stress model for copper JPreston et"al] J2003)). 
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Figure 10: Predicted values of yield stress from the Preston-Tonks-Wallace model. Please see the caption of Figure[5] 
for the sources of the experimental data. 
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Table 9: Comparison of the error in the yield stress predicted by the five flow stress models at various strain-rates and 
temperatures. 



1C111JJ. ^ I v ) 


Strain Rntp ( f<i\ 


Error 






7 A 


1V1 1 J V 


PTW (°/r>\ 




1 


IVlaX. 


^9 
9Z 




9 


z 


9 






\/Ti n 


A 


9 1 


1 n 


/I 


/; 
-o 






A/Iean 


0.2 


41 


.4 


0.2 


0.5 






IVlcUlall 


9 

J 


4.1 




fl 

w.u 


i i 
i . i 






ctH n™ 

OLU. UC\. 


f. 


7 




1 ^ 

1 .9 


9 7, 

Z.J 




8000 
oUUU 


IVlaX. 


1 1 

1 . 1 


J) 




1 9 
-1Z 


/; 
-u 






iviin . 


99 
-zz 


1 9 
-1Z 


91 
-Z 1 


9Q 
-zy 


9Q 






Mean 


-17 


-6 


-17 


-19 


-14 






\/T«H inn 
IvlCLllall 


-90 


.7 


-1 8 
1 o 


-1 8 

1 o 


-1 






OLU. 1VCV. 


u 




9 
z 


9 


A 


o / J 


9^00 


Max. 


7 
- / 


zLQ 


-9 


1 ^ 
13 


-9 






iviin . 


1 8 
-lo 


o 


9zl 


-9 


7 
- / 






Mean 


-13 


zo 


-19 


4 


zr 

-o 






Median 


-13 


25 


-16 


4 


-6 






Std. Dev. 


4 


16 


7 


7 


0.5 


1023 


1800 


Max. 


-16 


53 


3 


20 


-7 






Min. 


-30 


-3 


-22 


-7 


-13 






Mean 


-25 


17 


-13 


4 


-10 






Median 


-27 


11 


-17 


1.5 


-9 






Std. Dev. 


5 


21 


9 


10 


2 


1173 


0.066 


Max. 


93 


440 


149 


99 


7 






Min. 


39 


186 


119 


81 


3 






Mean 


64 


297 


132 


90 


5 






Median 


61 


275 


131 


92 


6 






Std. Dev. 


20 


93 


12 


6 


1.4 


1173 


960 


Max. 


-37 


50 


14 


24 


-13 






Min. 


-49 


-8 


-8 


-0.1 


-17 






Mean 


-45 


12 


-2 


9 


-15 






Median 


-47 


4 


-6 


6 


-14 






Std. Dev. 


4 


20 


8 


9 


1 



23 



Table 10: Comparison of the error in the yield stress predicted by the five flow stress models for a strain-rate of 4000/s. 
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Table 1 1 : Comparison of average "maximum" absolute (MA) errors in yield stresses predicted by the five flow stress 
models for various conditions. 
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model. The SCGL model again performs the worst. Finally, the low temperature tests (< 800 K) are predicted best 
by the PTW model. The other models also perform reasonably well under these conditions. 

From the above comparisons, the Preston-Tonks-Wallace and the Mechanical Threshold Stress models clearly stand 
out as reasonably accurate over the largest range of strain-rates and temperatures. To further improve our confidence 
in the above conclusions, we perform a similar set of comparisons with Taylor impact test data in the next section. 

Note that we could potentially recalibrate all the models to get a better fit to the experimental data and render the 
above comparisons void. However, it is likely that the average user of such models in computational codes will use 
parameters that are readily available in the literature with the implicit assumption is that published parameters 
provide the best possible fit to experimental data. Hence, exercises such as ours provide useful benchmarks for the 
comparative evaluation of various flow stress models. 



4 Taylor impact simulations 

The Taylor impact test dTavloJ (fl948)) was originally devised as a means of determining the dynamic yield strength 
of solids. The test involves the impact of a flat-nosed cylindrical projectile on a hard target at normal incidence. The 
test was originally devised to determine the yield strengths of materials a t high strain-rates. However, th at use of the 
test is limited to peak strains of around 0.6 at the center of the specimen Ijohnson and Holmauistl lll988h '). For higher 
strains and strain -rates, the Taylor test is more useful as a means of validating high strain-rate plasticity models in 
numerical codes dZerilli and Armstrong! i 19871) ). 

The attractiveness of the Taylor impact test arises because of the simplicity and inexpensiveness of the test. A 
flat-ended cylinder is fired on a target at a relatively high velocity and the final deformed shape is measured. The 
drawback of this test is that intermediate states of the cylinder are relatively difficult to measure. 

In this section, we compare the deformed profiles of Taylor cylinders from experiments with profiles that we obtain 
from our simulations. The experimental profiles are from the open literature and have been digitized at a high 
resolution. The errors in digitization are of the order of 2% to 5% depending on the clarity of the image. Our 
simulations use the Uintah code and the Material Point Method (seelAlandlEl. 

All our simulations are three-dimensional and model a quarter of the cylinder. We have used 8 material points per cell 
(64 material points per cell for simulations at 1235 K), a 8 point interpolation from material poin ts to grid, and a c ell 
spacing of 0.3 mm. A cell spacing of 0.15 mm gives essentially the same final deformed profile (lBanerieel ( E005bT) '). 
The anvil is modeled as a rigid material. Contact betwee n the cylinder and the anvil is assumed to be frictionless. The 
effect of frictional contact has been discussed elsewhere (Baneriee (2005b)). We have not included the effect of 
dama ge accum ulation due to void nucleation and growth in these simulations. Details of such effects can be found in 
Bane rieelj2005bl) . 

Our simulations were run for 150 /is - 200 /is depending on the problem. These times were sufficient for the 
cylinders to rebound from the anvil and to stop undergoing further plastic deformation. However, small elastic 
deformations continue to persist as the stress waves reflect from the surfaces of the cylinder. 

We have performed a systematic and extensive set of verification and validation tests to determine the accuracy of the 
Material Point Method and its implementation within Uintah ( iBaneried ll2004clJbll2005clalbh l. A number of 
materials and conditions have been explored in the process. We are, therefore, reasonably confident in the results of 
our simulations. 
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4.1 Metrics 



The systematic verification and validation of computational codes and the associated materi a l mod els requires the 
develo pment and utilization of appropriate comparison metrics (see Oberkam pf et a i1 d2002l) : lBabuska and Odenl 

In this section we discuss a few geometrical metrics that can be used in the context of Taylor impact tests. 
Other metrics such as the surface temperature and the time of impact may also be used if measured values are 
available. 

In most papers on the simulation of Taylor impact tests, a plot of the deformed configuration is superimposed on the 
experimental data and a visual judgement of accuracy is made. However, when the number of Taylor tests is large, it 
is not possible to present sectional/plan views for all the tests and numerical metrics are preferable. Some such 
metrics that have been used to compare Taylor impact tests are (see Figure ITT1 : 

1. The final length of the deformed cylinder (LA Jwilkins and Guinanldl973l) : lGustl ( ll982i) : IJones and Gillisl 
dl987l) : IJohnson and Holmauistl dl988l) : lHouse et alJdl995l) ). 

2. The di ameter of the mushroomed end of the cylinder (D j) djohnson and Holmauistl dl988l) : lHouse et all 
dl995l) V 

3. The length of the elastic zone in the cylinder (Xf) dJones and Gillisl dl987l) : lHouse et alJd 19951) 1. 

4. The bulge at a given distance from the deformed end (W/) djohnson and Holmauistl (fl988)). 

Contours of plastic strain have also been presented in a number of works on Taylor impact. However, such contours 
are not of much use when comparing simulations with experiments (though they are useful when comparing two 
stress update algorithms). 

The above metrics are inadequate when comparing the secondary bulges in two Taylor cylinders. We consider some 
additional geometrical metrics that act as a substitute for detailed pointwise geometrical comparisons between two 
Taylor test profiles. These are (see Figure ITTl : 

1. The final length of a axial line on the surface of the cylinder (L a f). 

2. The area of the cross-sectional profile of the deformed cylinder (Af). 

3. The volume of the deformed cylinder (Vf). 

4. The location of the centroid of the deformed cylinder in terms of a orthonormal basis with origin at deformed 
end (C xf , C y f). 

5. The moments of inertia of the cross section of the deformed cylinder about the basal plane (I x f) and an axial 
plane (I yf ). 

Higher order moments should also be computed so that we can dispense with arbitrary measures such as Wf. The 
numerical formulas used to compute the area, volume, centroid, and moments of inertia are given in AppendixO 

4.2 Experimental data 

In this section, we show plots of the experimentally determined values of some of the metrics discussed in the 
previous section. Quantities with subscript '0' represent initial values. The abscissa in each plot is a measure of the 
total energy density in the cylinder. The internal energy density has been added to the kinetic energy to separate the 
high temperature and low temperature data. 

Figure[21shows the length ratios (L f /Lq) for a number of Taylor impact tests. The figure indicates the following: 



26 



Table 12: Initial data for copper simulations. 



Test 


Material 


Initial 


Initial 


Initial 


Initial 


Source 






Length 


Diameter 


Velocity 


Temp. 








(L mm) 


(D mm) 


(V m/s) 


(To K) 




Cu-1 


OFHC Cu 


23.47 


7.62 


210 


298 


Wilkins and Guinan (1973) 


Cu-2 


ETPCu 


30 


6.00 


188 


718 


Gust ( 19821 


Cu-3 


ETPCu 


30 


6.00 


178 


1235 


Gust (1982) 



1 . The ratio (Lf/Lo) is essentially independent of the initial length and diameter of the cylinder. 

2. There is a linear relationship between the ratio (£// Lq) and the initial kinetic energy density. 

3. As temperature increases, the absolute value of the slope of this line increases. 

4. The deformation of OFHC (Oxygen Free High Conductivity) cannot be distinguished from that of ETP 
(Electrolytic Tough Pitch) copper from this plot. 

We have chosen to do detailed comparisons between experiment and simulation for the three tests marked with 
crosses on the figure. These tests represent situations in which fracture has not been observed in the cylinders and 
cover the range of temperatures of interest to us. 

The ratio of the diameter of the deformed end to the original diameter (D f /Do) for some of these tests is plotted as a 
function of the energy density in FigurefO] A linear relation similar to that for the length is observed. 

The volume of the cylinder should be preserved during the Taylor test if isochoric plasticity holds. Figure [T4l shows 
the ratio of the final volume to the initial volume (Vf /Vq) as a function of the energy density. We can see that the 
volume is preserved for three of the tests but not for the rest. This discrepancy may be due to errors in digitization of 
the profile. 

4.3 Evaluation of flow stress models 

In this section we present results from simulations of three Taylor tests on copper, compute validation metrics, and 
compare these metrics with experimental data. Tablell2lshows the initial dimensions, velocity, and temperature of the 
three specimens that we have simulated. All three specimens had been annealed before experimental testing. 

4.3.1 Test Cu-1 

Test Cu-1 is a room temperature test at an initial nominal strain-rate of around 9000/s. Figures [T5T a). (b), (c), (d), and 
(e) show the profiles computed by the JC, SCGL, ZA, MTS, and PTW models, respectively, for test Cu-1. 

The Johnson-Cook model gives the best match to the experimental data at this temperature (room temperature) if we 
consider the final length and the final mushroom diameter. All the other models underestimate the mushroom 
diameter but predict the final length quite accurately. The MTS model underestimates the final length. 

The time at which the cylinder loses all its kinetic energy (as predicted by the models) is shown in the energy plot of 
Figure [T5T f). The predicted times vary between 55 micro sees to 60 micro sees but are essentially the same for all the 
models. The total energy is conserved relatively well. The slight initial dissipation is the result of the artificial 
viscosity in the numerical algorithm that is used to damp out initial oscillations. 

In Figure^^we have seen that the final volume of the cylinder for test Cu-1 is around 5% larger than the initial 
volume. We assume that this error is due to errors in digitization. In that case, we have errors of + 1% for measures of 
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Figure 1 1 : Geometrical metrics used to compare profiles of Taylor impact specimens. 
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Figure 15: Computed versus experimental profiles for Ta^Jpr test Cu-1 and the computed energy versus time profile. 
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Figure 16: Comparison of error metrics for the five models for Taylor test Cu-1. 



length and errors of +2% for measures of area in the experimental profile. Moments of inertia of areas are expected to 
have errors of around 7%. 

The error metrics for test Cu-1 are shown in Figure^] The final length (Lf) is predicted to within 3% of the 
experimental value by all the models. The Johnson-Cook and Preston-Tonks-Wallace models show the least error. 

The length of the deformed surface of the cylinder (L a f) is predicted best by the Johnson-Cook and 
Steinberg-Cochran-Guinan-Lund models. The other models underestimate the length by more than 5%. The final 
mushroom diameter (Df) is underestimated by 5% to 15%. The Johnson-Cook model does the best for this metric, 
followed by the Mechanical Threshold Stress model. The width of the bulge (Wf) is underestimated by the 
Johnson-Cook and SCGL models and accurately predicted by the ZA, MTS, and PTW models. The length of the 
elastic zone (Xf) is predicted to be zero by the SCGL, ZA, MTS, and PTW models while the Johnson-Cook model 
predicts a value of 1.5 mm. Moreover, an accurate estimate of Xf cannot be made from the experimental profile for 
test Cu-1. Therefore we do not consider this metric of utility in our comparisons for this test. 

From Figure^]we see that the predicted area of the profile (Af) is within 3% of the experimental value for all the 
models. The SCGL model shows the least error in this metric. If we decrease the experimental area by 2% (in 
accordance with the assumed error in digitization), the Johnson-Cook and PTW models show the least error in this 
metric. 

The predicted final volume of the cylinder is around 0.8% larger than the initial volume showing that volume is not 
preserved accurately by our stress update algorithm. The error in digitization is around 5%. That gives us a uniform 
error of 5% between the experimental and computed volume (Vf) as can be seen in Figure[T6l 

The locations of the centroids (C x f, C y f) provide further geometric information about the shapes of the profiles. 
These are the first order moments of the area. The computed values are within 2% of experiment except for the MTS 
model which shows errors of -4% for C x f and +6% for C y f. 

The second moments of the area are shown as I x f and I y f in Figure^] The error in I x f tracks and accentuates the 
error in Lf while the error in I y f tracks the error in Df . The width of the bulge is included in this metric and it can 
be used the replace metrics such as Lf, Df, and Wf for the purpose of comparison. We notice this tracking behavior 
when the overall errors are small but not otherwise. 

We have also plotted the arithmetic mean of the absolute value of the errors in each of the metrics to get an idea about 
which model performs best. The average error is the least (2.5%) for the Johnson-Cook model, followed by the ZA 
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and PTW models (3.5%). The MTS model shows an average error of 4% while the SCGL model shows the largest 
error (5%). If we subtract the digitization error from the experimental values, these errors decrease and lie in the 
range of 2% to 3%. 

In summary, all the models predict profiles that are within the range of ex perimental variation for the test at room 
temperature. Additional simulations at higher strain-rates ( iBanerieel I l2005bl) ) have confirmed that all the models do 
well for room temperature simulations for strain rates ranging from 500 /s to 8000 /s. We suggest that the simplest 
model should be used for such room temperature simulations and our recommendation is the Zerilli-Armstrong 
model for copper. 

4.3.2 TestCu-2 

Test Cu-2 is at a temperature of 718 K and the initial nominal strain-rate is around 6200/s. Figures [T7F a). (b), (c), (d), 
and (e) show the profiles computed by the JC, SCGL, ZA, MTS, and PTW models, respectively, for test Cu-2. 

In this case, the Johnson-Cook model predicts the final length well but overestimates the mushroom diameter. The 
SCGL model overestimates the length but predicts the mushroom diameter well. The ZA model predicts the overall 
profile remarkably well except for the mushroom diameter. The MTS model slightly overestimates both the final 
length and the mushroom diameter. The PTW model also performs similarly, except that the error is slightly larger 
than that for the MTS model. 

The energy plot for test Cu-2 is shown in Figure fTTl f). In this case, the time of impact predicted by the JC and ZA 
models is around 100 micro sees while that predicted by the SCGL, MTS, and PTW models is around 90 micro sees. 

The error metrics for test Cu-2 are shown in Figure^] In Figure[2]we have seen that the deformed volume 
computed from the digitized profile is almost exactly equal to the initial volume for test Cu-2. The digitization error 
can be neglected in this case. 

The least error in the predicted final length (L /) is for the ZA model followed by the JC model. The SCGL model 
shows the largest error in this metric (7%). The MTS and PTW models overestimate the final length by around 6%. 
The value of L a f is predicted to within 2% of the experimental value by the ZA model. The corresponding errors in 
the other models vary from 6% (MTS) to 9% (SCGL). The mushroom diameter is overestimated by all models. The 
JC model overestimates this metric by more than 30%. The ZA and PTW models overestimate Df by 17% to 19%. 
The MTS model overestimates Df by 12%. The SCGL model does best with an error of 7%. The width of the bulge 
is underestimated by all the models with errors varying between 5% (ZA) tp 9% (JC). 

The final area (^4/) is predicted almost exactly by the JC model. The ZA model underestimates the area by 1% while 
the errors in the other models vary from 2% to 4%. The error in the final volume is less than 1% for all the models. 

The location of the centroid is predicted best by the Johnson-Cook model followed by the ZA model. Both the MTS 
and PTW models underestimate C x f by 2% and overestimate C y f by 5%. The SCGL model shows the largest error 
for this metric. 

For the second order moments I x f, the smallest error is for the Johnson-Cook model followed by the ZA model. The 
largest errors are from the SCGL model. The MTS and PTW models overestimate this metric my 15%. The PTW 
model predicts I y f the best, followed by the MTS model showing that the overall shape of the profile is best 
predicted by these models. The Johnson-Cook and SCGL models show the largest errors in this metric. 

On average, the ZA model performs best for test Cu-2 at 718 K with an average error of 4%. The MTS model shows 
an average error of 5% while the JC and PTW models show errors of approximately 6%. The SCGL model, with an 
average error of approximately 7%, does the worst. 
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Figure 17: Computed and experimental profiles for Taylor test Cu-2 and the computed energy-time profile. 




Figure 18: Comparison of error metrics for the five models for Taylor test Cu-2. 



4.3.3 Test Cu-3 

Test Cu-3 was conducted at 1235 K and at a initial nominal strain-rate of approximately 6000 /s. Figures [T9F a). (b), 
(c), (d), and (e) show the profiles computed using the JC, SCGL, ZA, MTS, and PTW models, respectively. The 
complete experimental profile of the cylinder was not available for this test. 

The Johnson-Cook model fails to predict a the deformation of the cylinder at this temperature and the material 
appears to flow along the plane of impact. The SCGL model predicts a reasonably close value of the final length. 
However, the low strain-rate part of the SCGL model behaves in an unstable manner at some levels of discretization 
for this test and should ideally be discarded in high strain-rate simulations. The ZA model overestimates the final 
length as does the MTS model. The PTW model predicts a final length that is closer to experiment but does not show 
the bulge that is characteristic of hardening. This can be seen from the tendency of the model to saturate prematurely 
as discussed in the section on one-dimensional tests. 

The energy plot for test Cu-3 is shown in Figure [T9F f). For this test, the JC model predicts a time of impact greater 
than 250 micro sees while the rest of the models predict values between 120 micro sees and 130 micro sees. The 
reason for the anomalous behavior of the JC model is that the rate dependence of the yield stress at high temperature 
is severely underestimated by the JC model. The nominal strain-rate is around 5000/s for this test at which the yield 
stress should be considerably higher than the 50 MPa that is computed by the JC model. 

We do not have the final profile of the sample for this test and hence cannot compare any metrics other than the final 
length. The final length is predicted most accurately by the SCGL model with an error of 10%, followed by the PTW 
model (error 15%) and the ZA model (error 20%). The Johnson-Cook model shown an error of more than 90%. 

These three sets of tests show that the performance of the models deteriorates with increasing temperature. However, 
on average all the models predict reasonably accurate profiles for the Taylor impact tests. The choice of the model 
should therefore be dictated by the required computational efficiency and the conditions expected during simulations. 



5 Summary and conclusions 

We have compared five flow stress models that are suitable for use in high strain-rate and high temperature 
simulations using one-dimensional tension/compression tests and Taylor impact tests. We have also evaluated the 
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Figure 19: Computed and experimental profiles for Taylor test Cu-3 and the computed energy-time profiles. 



associated models for shear modulus, melting temperature, and the equation of state. We observe that during the 
simulation of large plastic deformations at high strain-rates and high temperatures, the following should be taken into 
consideration: 



1. The specific heat can be assumed constant when the range of temperatures is small. However, at temperatures 
below 250 K or above 750 K, the room temperature value of the specific heat may not be appropriate. 

2. The Mie-Griineisen equation of state that we have used is valid only up to compressions of 1.3. A higher order 
approximation should be used if extreme pressures are expected during the simulation. We note that care 
should be exercised when this equation of state is used for states of large hydrostatic tension. 

3. The physically-based Burakovsky-Preston-Silbar melt temperature model should be used when the material is 
not well characterized. However, the Steinberg-Cochran-Guinan model is quite accurate for copper and slightly 
less computationally expensive. 

4. The Nadal-Le Poac shear modulus model is the most robust model for high-temperature high-pressure 
applications. However, this model (as well as the Steinberg-Cochran-Guinan model) can be unreliable at high 
hydrostatic tensions. The MTS shear modulus model is the best at temperatures less than 200 K though it is 
inaccurate at high pressures. 

5. The Johnson-Cook (JC) flow stress model predicts that the strain-rate dependence of the yield stress is 
insignificant at high temperatures. This is inaccurate and the model should be recalibrated for high temperature 
applications. 

6. The contribution of the rate-dependent part of the Steinberg-Cochran-Guinan-Lund (SCGL) model to the total 
flow stress is small enough to be insignificant because it is limited by a constant Peierls stress. The model 
makes the assumption that da y /dT = d/i/dT. This assumption clearly does not hold for copper. If the model 
is corrected, the high strain-rate yield stresses are predicted quite accurately by the model. 

7. The Zerilli-Armstrong (ZA) model is quite accurate at low strain-rates. The temperature-dependence of the 
yield stress predicted by the model is larger than the actual value at high strain-rates with some exceptions. 

8. The Mechanical Threshold Stress (MTS) model shows the correct temperature dependence of the yield stress 
and high temperatures and strain-rates. The model is also quite accurate at low strain-rates. However, the 
stage-IV hardening at large tensile strains is not predicted by the model. 

9. The Preston-Tonks-Wallace (PTW) model is also quite accurate at both low and high strain-rates and at high 
temperatures. However, the predicted yield stress saturates at a tensile strains of 0.4 and the model represents 
stage-IV hardening worse than the Mechanical Threshold Stress model. 

10. The overall error for one-dimensional tests is least for the PTW model. For the tension tests, the MTS model is 
the most accurate while the PTW model is the most accurate for compression tests. The MTS model is the 
most accurate at high strain-rates while the PTW model is the most accurate at low strain-rates. The PTW 
model is the most accurate for both high- and low-temperature tests at all strain-rates. 

1 1 . For comparing Taylor impact tests, geometric moments of the profile are better measures of the shape than 
arbitrary measures such as the width of the bulge at a fixed distance from the end. 

12. At room temperature, all the flow stress models predict the final profile of a Taylor impact cylinder quite 
accurately. At 718 K, the ZA model is the most accurate, followed by the MTS model. At 1235 K, the JC 
model fails to predict the final profile while the modified SCGL and the PTW models predict the final length 
most accurately. The PTW model fails to predict the bulge due to hardening at this temperature due to early 
saturation of the yield stress. 



We conclude that no single model can accurately predict the yield stress for the full range of conditions considered in 
this study. Care should be exercised when choosing the models to be used in a particular simulation. If needed, the 
model should be recalibrated for the range of conditions of interest. 
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A The Material Point Method 

The Material Point Method (MPM) Sulsk v et alJ ( 1 1994 1 1995b is a particle method for solid mechanics simulations. 
In this method, the state variables of the material are described on Lagrangian particles or "material points". In 
addition, a regular, structured Eulerian grid is used as a computational scratch pad to compute spatial gradients and to 
solve the governing conservation equations. The explicit time-stepping version of the Material Point Method is 
summarized below. 

It is assumed that an particle state at the beginning of a time step is known. The mass (m), external force (/ ext ), and 
velocity (v) of the particles are interpolated to the grid using the relations 

m g = y ^2,S gp m p , v g = {l/m g )^2s gp m p v p , f g xt = ^2s gp f p xt (36) 
p p p 

where the subscript 'g' indicates a quantity at a grid node and a subscript 'p' indicates a quantity on a particle. The 
symbol indicates a summation over all particles. The quantity (S gp ) is the interpolation function of node g 
evaluated at the position of particle p. 

Next, the velocity gradient at each particle is computed using the grid velocities using the relation 

Vv p = G gp v g (37) 

a 

where G gp is the gradient of the shape function of node g evaluated at the position of particle p. The velocity 
gradient at each particle is used to determine the Cauchy stress (cr p ) at the particle using a stress update algorithm. 

The internal force at the grid nodes (/"') is calculated from the divergence of the stress using 

ff = Y J G gp <r p V p (38) 

p 

where V p is the particle volume. 

The equation for the conservation of linear momentum is next solved on the grid. This equation can be cast in the 
form 

m g a g = ff - ff (39) 
where a g is the acceleration vector at grid node g. 

The velocity vector at node g is updated using an explicit (forward Euler) time integration, and the particle velocity 
and position are then updated using grid quantities. The relevant equations are 

v g {t + At) =v g (t)+a g At (40) 

v p (t + At) = v p (t) +^2s gp a g At ; x p (t + At) = x p (t) +^ S gp v g At (41) 

9 9 

The above sequence of steps is repeated for each time step. The above algorithm leads to particularly simple 
mechanisms for ha ndling contact. Deta ils of the interpolants a nd the contact algorithms can be found in 
iBardenhagen et alJ J2001I) andlBardenhagen and Kober<2004). The implicit version of the Material Point method can 



be found in lGuilkev and Weis 
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B Stress Update Algorithm 



A modified form of a hypoelastic - plastic, semi-imp l icit elastic-plastic stress update algo r ithm dNemat-NasseJ Jl99ll) : 
iNemat-Nasser and Chungl Jl992l) : [wang and Atiuril Jl994 : lMaudlin and Schiferll ( fl996l) : IZocher et alJ l l200ol) ) has 
been used for the stress update in the simulations presen ted in this paper. D etails of the socket-based objected 
oriented form of the complete algorithm can be found in lBaneried l2005b). We provide a brief description here for 
the sake of completeness. 

Following iMaudlin and Schiferll (H 996). the rotated spatial rate of deformation tensor (d) is decomposed into an 
elastic part (d e ) and a plastic part (d p ) 

d = d e + d p (42) 
If we assume plastic incompressibility ( tr(d p ) = 0), we get 

r] = V e + V P (43) 

where r), rf, and rf are the deviatoric parts of d, d e , and d p , respectively. For isotropic materials, the hypoelastic 
constitutive equation for deviatoric stress is 

s = 2n(rj - rf) (44) 

where s is the deviatoric part of the stress tensor and /1 is the shear modulus. We assume that the flow stress obeys the 
Huber-von Mises yield condition 



/ := y 2 11*11 - °v < or ' F:=-s:s-ai<0 (45) 
where a y is the flow stress. Assuming an associated flow rule, and noting that d v = r/ p , we have 

rf = d p = = = 3As (46) 

OCT OCT 

where er is the stress. Let u be a tensor proportional to the plastic straining direction, and define 7 as 

u = v / 3 7 A r ; j:=VSA\\s\\ =^> - ( u = 3As (47) 

ll s ll 

Therefore, we have 

rf = ju; s = 2/j,(r) — ju) (48) 
From the consistency condition, if we assume that the deviatoric stress remains constant over a timestep, we get 

s : n 

7 = ' (49) 



which provides an initial estimate of the plastic strain-rate. To obtain a semi-implicit update of the stress using 
equation d48i . we define 

r 2 := : a = o* (50) 
Taking a time derivative of equation i5Q\ gives us 

V2f = V3^-4 (51) 

HI 

Plugging equation J5 1 1 into equation d48> 9 we get 

f = V2~n(u : r) — 7M : u) = v2;u(d — 37) (52) 
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where d = u : rj. If the initial estimate of the plastic strain-rate is that all of the deviatoric strain-rate is plastic, then 
we get an approximation to 7, and the corresponding error (7^) given by 

d d 



(53) 



The incremental form of the above equation is 

A7 = - A 7er (54) 
Integrating equation d52i from time t n to time t n+ \ = t n + At, and using equation i54\ we get 

Tn+i = r n + V2fi(d*At - 3A 7 ) = T „ + 3V2^A lei (55) 
where d* is the average value of d over the timestep. Solving for A7 el gives 

. T n+ l - T n V2(Jy ~ V3\\s n \\ - 

A%r = ^v%r= — ^ — (56) 

The direction of the total strain-rate (u v ) and the direction of the plastic strain-rate (tt s ) are given by 

\\v\\ \\s\\ 

Let 9 be the fraction of the time increment that sees elastic straining. Then 

d* - 37„ 



d* 

where j n = d n /3 is the value of 7 at the beginning of the timestep. We also assume that 



(58) 



d* = V3»7 : [(1 - 9)u r ' + -{u r < + u s )} (59) 
Plugging equation d58l into equation d59t we get a quadratic equation that can be solved for d* as follows 

-^=(d*) 2 ~(v-u s + IMIK + 37n(»7 : ~ Nl) = (60) 

The real positive root of the above quadratic equation is taken as the estimate for d. The value of A7 can now be 
calculated using equations J54> and i56\ . A semi-implicit estimate of the deviatoric stress can be obtained at this 
stage by integrating equation ( 1481 9 

s n +i = s n + I v At - V3A 7 /^ ) (61) 
= s n + 2fi (r)At - A A7 ^±i^ (62) 

Solving for S„+i, we get 

c trial 

s n +l = -r- (63) 

1 + 3y/2n— 

<7y 

where s t "f 1 = s n + 2/.iAtr]. A final radial return adjustment is used to move the stress to the yield surface 

[2 Sn+l 

S " +1 = V3^P^ (64) 

A pathological situation arises if j n = u n : rj n is less than or equal to zero or A7 er > 4- At. This can occur is the 
rate of plastic deformation is small compared to the rate of elastic deformation or if the timestep size is too small 
(see Nemat ^Nasser and Chund(ll992l)'). In such situ ations, we use a locally implicit stress update that uses Newton 
iterations (as discussed in lsrmoandHughesl ( ll998l) . page 124) to compute s. 
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C Computation of Metrics 



The length of the elastic zone after deformation (Xf) is determined by checking the deformed diameter with the 
original diameter of the cylinder. If the difference is greater than 0.003 mm, plastic deformation is assumed to have 
taken place. The value of Xf is the distance from the free end of the cylinder to the first point from the free end 
where the above criterion is met. 

Let the closed polygon representing the final profile of the Taylor cylinder be given by 

P = pi,P2,P3, ■■■■,Pn,Pn+i = Pi, where n is the number of vertices of the polygon. We assume that the points are 
ordered in the counter-clockwise direction. Each point pi has a pair of coordinates (xu yi). 



Then, the area of the profile (Af) is given by 



1 

A f = o Vl+1 ~ Xi+1 yA 

1 i=i 



(65) 



The centroid of the profile is given by 



l ,A 

C x f = y fa yi+i - Xi+i yi){xi + Xi+i) (66) 



6A ftt 
1 " 

Cyf = g^X^ 1 Vl+1 - Vi){Vi + Vi+i) • 

■* i=l 



(67) 



The volume of the deformed cylinder is given by the Pappus theorem. The formula for the volume is 

V f = 2nC xf A f . (68) 

The moments of inertia are computed by converting the volume integral into a surface integral over the boundary of 
the profile. The resulting formulas for the moments of inertia are 

1 " 

T xf = -j^2( x *+i - x i)(yi+i + y l )(vhi + vf) ( 69 ) 
i " 

hi = ~[j^(y^ - Vi){xi+i + Xi)(Xi +1 + xf) . (70) 
i=i 
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